%% process_error_NS_s
function [error_u, error_p] = process_error_NS_s(mesh, fespace1, fespace2, pde, U, error_type)
Nb1 = fespace1.Nb;
switch error_type
    case "Linf"
        err_u1 = 0;
        err_u2 = 0;
        err_p = 0;
        for n = 1:mesh.Ne
            V = mesh.node(mesh.elem(n,:),:);
            [int_points_loc, ~] = integral_loc(V, "2d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "x");
            psi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace2.basis_type, "x");
            u1h = phi * U(fespace1.elem(n,:));
            u2h = phi * U(fespace1.elem(n,:)+Nb1);
            ph = psi * U(fespace2.elem(n,:)+2*Nb1);
            f1 = max(abs(pde.u1(int_points_loc(:,1),int_points_loc(:,2)) - u1h));
            f2 = max(abs(pde.u2(int_points_loc(:,1),int_points_loc(:,2)) - u2h));
            f3 = max(abs(pde.p(int_points_loc(:,1),int_points_loc(:,2)) - ph));
            if f1 > err_u1
                err_u1 = f1;
            end
            if f2 > err_u2
                err_u2 = f2;
            end
            if f3 > err_p
                err_p = f3;
            end
        end
        error_u = max(err_u1, err_u2);
        error_p = err_p;
    case "L2"
        err_u1 = 0;
        err_u2 = 0;
        err_p = 0;
        for n = 1:mesh.Ne
            V = mesh.node(mesh.elem(n,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(V, "2d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "x");
            psi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace2.basis_type, "x");
            u1h = phi * U(fespace1.elem(n,:));
            u2h = phi * U(fespace1.elem(n,:)+Nb1);
            ph = psi * U(fespace2.elem(n,:)+2*Nb1);
            f1 = (pde.u1(int_points_loc(:,1),int_points_loc(:,2)) - u1h).^2;
            f2 = (pde.u2(int_points_loc(:,1),int_points_loc(:,2)) - u2h).^2;
            f3 = (pde.p(int_points_loc(:,1),int_points_loc(:,2)) - ph).^2;
            r1 = int_weights_loc * f1;
            r2 = int_weights_loc * f2;
            r3 = int_weights_loc * f3;
            err_u1 = err_u1 + r1;
            err_u2 = err_u2 + r2;
            err_p = err_p + r3;
        end
        error_u = sqrt(err_u1 + err_u2);
        error_p = sqrt(err_p);
    case "H1"
        err_u1x = 0;
        err_u1y = 0;
        err_u2x = 0;
        err_u2y = 0;
        err_px = 0;
        err_py = 0;
        for n = 1:mesh.Ne
            V = mesh.node(mesh.elem(n,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(V, "2d");
            phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "dx");
            phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "dy");
            psix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace2.basis_type, "dx");
            psiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace2.basis_type, "dy");
            u1hx = phix * U(fespace1.elem(n,:));
            u1hy = phiy * U(fespace1.elem(n,:));
            u2hx = phix * U(fespace1.elem(n,:)+Nb1);
            u2hy = phiy * U(fespace1.elem(n,:)+Nb1);
            phx = psix * U(fespace2.elem(n,:)+2*Nb1);
            phy = psiy * U(fespace2.elem(n,:)+2*Nb1);
            f1x = (pde.u1dx(int_points_loc(:,1),int_points_loc(:,2)) - u1hx).^2;
            f1y = (pde.u1dy(int_points_loc(:,1),int_points_loc(:,2)) - u1hy).^2;
            f2x = (pde.u2dx(int_points_loc(:,1),int_points_loc(:,2)) - u2hx).^2;
            f2y = (pde.u2dy(int_points_loc(:,1),int_points_loc(:,2)) - u2hy).^2;
            f3x = (pde.pdx(int_points_loc(:,1),int_points_loc(:,2)) - phx).^2;
            f3y = (pde.pdy(int_points_loc(:,1),int_points_loc(:,2)) - phy).^2;
            r1x = int_weights_loc * f1x;
            r1y = int_weights_loc * f1y;
            r2x = int_weights_loc * f2x;
            r2y = int_weights_loc * f2y;
            r3x = int_weights_loc * f3x;
            r3y = int_weights_loc * f3y;
            err_u1x = err_u1x + r1x;
            err_u1y = err_u2y + r1y;
            err_u2x = err_u2x + r2x;
            err_u2y = err_u2y + r2y;
            err_px = err_px + r3x;
            err_py = err_py + r3y;
        end
        error_u = sqrt(err_u1x + err_u1y + err_u2x + err_u2y);
        error_p = sqrt(err_px + err_py);
    otherwise
        error("Invalid error type.");
end
end